Reaching large lengths and long times in polymer dynamics simulations 



A. van Heukelum and G. T. Barkema 
Institute for Theoretical Physics, Utrecht University, 
Leuvenlaan 4, 3584 CE Utrecht, The Netherlands 

A lattice model is presented for the simulation of dynamics in polymeric systems. Each polymer 
is represented as a chain of monomers, residing on a sequence of nearest-neighbor sites of a face- 
centered-cubic lattice. The polymers are self- and mutually avoiding walks: no lattice site is visited 
by more than one polymer, nor revisited by the same polymer after leaving it. The dynamics 
occurs through single- monomer displacements over one lattice spacing. To demonstrate the high 
computational efficiency of the model, we simulate a dense binary polymer mixture with repelling 
nearest-neighbor interactions between the two types of polymers, and observe the phase separation 
over a long period of time. The simulations consist of a total of 46, 080 polymers, 100 monomers 
each, on a lattice with 13,824,000 sites, and an interaction strength of O.lfcsT. In the final two 
decades of time, the domain-growth is found to be d(t) ~ t ly/3 , as expected for a so-called "Model B" 
system. 



I. INTRODUCTION 

Lattice models are widely used as a theoretic tool to 
study polymer solutions. Most notably, Flory 0, and 
Huggins 0, introduced a mean- field theory of a partic- 
ular simple model which has proven powerful enough to 
explain solubility properties and liquid-liquid phase sep- 
aration of polymer mixtures. An extension to the model, 
describing the dynamics of such systems, has been intro- 
duced by Cahn and Hilliard HGf and Hillert 0. Some 
of the earliest computer simulations of lattice polymers 
were performed in 1962 by Verdier and Stockmayer 
who studied the dynamics of polymers, modeled by self- 
avoiding walks on square and cubic lattices. Kolinski et 
al. 0, added two-monomer moves to the dynamics 
and simulated a homopolymeric melt. 

To describe the polymer dynamics in a melt with more 
detail than the random-walk approach, Carmesin intro- 
duced in 1988 the bond fluctuation model ^ n the 
two-dimensional version of this model, each monomer 
occupies four (2 x 2) lattice sites of a square lattice. 
Multiple occupation of lattice sites is not permitted. 
Monomers adjacent in the chain are connected by bonds 
with lengths between 2 and The dynamics of 
the polymer chain consists of the displacement of sin- 
gle monomers to nearest-neighbor lattice sites, restricted 
by the constraints on bond length and excluded volume. 
The bond-fluctuation model was extended to three di- 
mensions m 1991, by Deutsch and Binde r!^ . This 
model has been used by several groups 0, Il4j| to simu- 
late dense polymer melts. 

In 1971, de Gennes proposed that the main mech- 
anism of polymer dynamics in gels and dense polymer 
mixtures is reptation, i.e., diffusion of stored length. To 
verify the theoretical predictions of de Gennes, Evans and 
Edwards [lfl ] introduced the cage model, in this model, a 
single polymer is simulated as a random walk on a square 
or cubic lattice. In the limit of tight gels, the dynamics 
reduces to the re-orientation of pairs of connected seg- 
ments pointing in opposite direction ("kinks"). The cage 
model has been used for the study of diffusion and relax- 



ation times of a single polymer in a gei miniiiii3, 

as well as for star polymers |2(j • Extending the model to 
include cxcludcd-volume effects or interactions between 
polymers is not easy. 



An alternative model for a single reptating polymer 
was proposed by Rubinstein 21]. In this model, known 
as the repton model, the polymer is represented as a chain 
of monomers, residing on the sites of a square or cubic 
lattice. Monomers adjacent along the chain reside in ci- 
ther nearest-neighbor lattice sites, or in the same site. 
The latter case corresponds to the presence of stored 
length, a central concept in de Gennes' description of en- 
tangled polymer dynamics. The dynamics of the polymer 
is strictly limited to reptation: the only kind of moves al- 
lowed to a monomer in the interior of the chain is that it 
can move to an adjacent lattice site, provided that one of 
its adjacent neighbors is already on the site to which it 
is moving, and the other is in the site it leaves. Per unit 
of time, each monomer attempts to move in each direc- 
tion statistically once. The repton model is illustrated in 
Figure^ The re pto n model has been studied extensively 
by one of us JHHl. 



The repton model can be extended to simulate poly- 
mers with self-avoiding walk statistics, as well as many 
interacting polymers; this is the route followed in this 
manuscript. In the next section, we discuss the exten- 
sions to the repton model which are required to describe 
equilibrium and dynamical properties of dense polymer 
solutions and melts. Since the dynamics of polymeric sys- 
tems occurs on long time-scales and large length-scales, 
computational efficiency is of paramount importance. We 
therefore discuss in detail how the resulting extended rep- 
ton model can be implemented with high computational 
efficiency, exploiting a technique known as multispin cod- 
ing. At the end of this manuscript, we show the useful- 
ness of our work by presenting simulation results on the 
phase separation of a dense binary polymer mixture. 
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FIG. 1: The repton model on a square lattice. The monomers 
in the polymer chain are connected by bonds of zero or unit 
length and each monomer is assigned to a face of the lattice, 
such that adjacent monomers in the chain occupy either the 
same or nearest-neighbor lattice sites. If a monomer is con- 
nected to one neighbor by a bond with zero length, and to 
its other neighbor by a bond with unit length, it can join its 
other neighbor. By doing so, the bond with zero length — a 
representation of "stored length" — diffuses along the chain. 
For instance, in the configuration depicted here, monomer 5 
can join its neighbor 6. The ends of the polymer allow fluc- 
tuations of the total amount of stored length. 




FIG. 2: Extended repton model in two dimensions. Sideways 
motion of a single monomer, preserving the constraint that 
bonds are either zero-length or unit-length, is only possible 
on lattices that contain loops of three sites, i.e., triangular 
(FCC) instead of the more usual square (cubic) lattice. The 
excluded volume constraint is implemented by allowing only 
adjacent monomers in a polymer to occupy the same lattice 
site. Phase separation is induced by nearest-neighbor lattice 
site interactions. 



II. EXTENSIONS OF THE REPTON MODEL 
TO SIMULATE POLYMER MELTS 

The repton model as proposed by Rubinstein simulates 
a single polymer obeying random-walk statistics, with 
dynamics limited to the diffusion of stored length. Two 
extensions to the repton model arc needed to use it for the 
simulation of entangled chains: the chains have to obey 
volume exclusion and on the other hand, the constraint 
that disallows sideways motion has to be lifted. 

The repton model describes "stored length" as two or 
more monomers on the same lattice site. Therefore, we 
cannot limit the polymers to self-avoiding-walk statistics 
by disallowing two monomers to occupy the same lattice 
site, without destroying the reptation dynamics that so 
successfully described polymer-diffusion in a gel. As rep- 
tation is also expected to play an important role in the 
melt, we should keep this kind of dynamics. The solu- 
tion to this problem is to limit only the contours of the 
polymers to self-avoiding-walk statistics. Thus, multi- 
ple occupation of a lattice site is allowed only for two or 
more monomers which are adjacent along the chain. One 
convenient side-effect of this choice is that since the rep- 
tation moves in the interior of the chain do not affect the 
tube, they never cause a violation of volume exclusion. 

Sideways motion is implemented by also allowing 
monomers to move by a single lattice spacing if this 
changes the tube of the polymer. Single-monomer moves 
work well on lattice structures that contain loops of three 
sites, like the triangular and face-centered-cubic (FCC) 



lattice; on lattices without such three-site loops, like the 
square and cubic lattices, the sideways movement can 
lead to lots of unwanted artifacts or, alternatively, must 
be implemented by allowing bonds to extend to next- 
nearest neighbors. 

The repton model is often studied in its projected ver- 
sion, in which the state of the polymer is characterized 
by the bonds Sj, i = 1 . . . N — 1, between the monomers 
i = 1 . . . N, where Si = Xi+i — Xi is equal to zero or ±1. 
Usually all allowed moves are attempted with the same 
(unit) rate, with the consequence that the three values 
for Si occur with equal probability. The density of stored 
length (probability that Si is zero) is then 1/3. To obtain 
the same density of stored length in higher-dimensional 
repton models (before projection), the rates for moves 
which decrease the stored length should equal 2/z times 
those of the reverse moves, in which z is the lattice co- 
ordination number: 4, 6 and 12 for the square, cubic or 
triangular, and FCC lattices, respectively. To conserve 
equilibrium properties, sideways moves which decrease 
the density of stored length should also be attempted 
with rates which are 2/z times the reverse moves. 

Figurc|21shows model-polymers on a triangular lattice. 
In the upper polymer in the figure, reptation moves are 
possible for monomers 2, 4, 6, 9, 10, and 11; sideways 
moves can be made by monomers 2, 4, 6, 8, 9, 10, and 11. 
Monomers 3, 5, and 7 cannot move. The end-monomers 
1 and 12 can move to any empty nearest-neighbor site. 
In the lower polymer, reptation moves are possible for 
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monomers 3, 5, 6, 10, and 11; sideways moves can be 
made by monomers 3, 4, 5, 6, 7, 10, and 11. Monomers 
2, 8, and 9 cannot move. The end-monomer 1 can move 
to any empty nearest-neighbor site, while end-monomer 
12 can only move to the site occupied by monomer 11. 

III. MULTISPIN CODING 

Multispin coding is a programming technique which 
makes use of the low-level bit manipulation functions 
provided by the CPU of our computer, effectively do- 
ing many simple computations in parallel. These func- 
tions include logical operations like AND (A), OR (V) 
exclusive-or (XOR, ffi), negation (NOT, ->), as well as 
bitwise shifts (<C, integer multiplication by a power of 2; 
3>, integer division by a power of 2); and arithmetic op- 
erations like addition and multiplication. For a general 
introduction to multispin coding we refer to Ref. p^j . 
Chapter 15. Here, we will explain this technique anal- 
ogous to Ref. 2J5j, as applied to the projected repton 
model. Next we will discuss the implementation of this 
technique in our model for polymer melts. 

A. repton model 

In the projected repton model, the state of the chain 
with monomers i = 1 . . . N is specified by the set of their 
coordinates {x\, . . . , xn}- Monomers which are adjacent 
along the chain have to reside in either the same lattice 
site, or nearest neighbors, with the consequence that the 
step Si = Xi+i — Xi can only take the values ±1 or 0. An 
alternative way to store the polymer state is therefore to 
store the coordinate x\ of the first monomer, and to store 
for step i = 1 ... N — 1 two arrays of bits (Ai, Bi) that 
can take the combination (0,0) if s$ = 0, (0, 1) if s$ = 1, 
and (1, 0) if S{ = —1. 

In one unit of time, each monomer attempts to move 
in each direction statistically once. Thus, 2N elemen- 
tary moves, i.e. one monomer attempting to move in one 
direction, should be attempted per unit of time. 

It is a property of the repton model, that if monomer 
i in the interior of the chain can join its neighbor i + 1, 
it cannot join its neighbor i — 1, and vice versa. We 
exploit this by attempting these two elementary moves 
simultaneously. The first and last monomer can move 
to both directions, if they are located in the same site 
as their neighbor along the chain, so this trick does not 
apply for these two monomers; they are therefore selected 
with a double probability. 

In detail, each monomer in the interior of the chain 
is selected with a probability l/(iV + 2). If the selected 
monomer i is located in the same site as exactly one of its 
two adjacent neighbors, it will join the other neighbor. 
The bits that contain information on whether a move is 
possible or not, are {Ai, Bi) and (Aj, Bj) with j = i — 1. 
A move is possible if one of the two pairs is (0, 0) and 
the other is not. If the move is carried out, (Ai, Bi) and 



(Aj, Bj) are exchanged. A sequence of logical operations 
that achieves this is: 

y = (AiV Bi)®(AjV Bj) 

m A = (Ai ® Aj) Ay 

m B = (Bi ® Bj) Ay 

A\ = Ai®m A (1) 

B[ = Bi®m B 

A'- = Aj © m A 

B'j = Bj © m B 

The first and last monomers are selected, each, with 
a probability 2/(N + 2); twice the probability to select 
a specific monomer in the interior of the chain. After 
the selection of an end monomer, the intended direction 
is also randomly selected, with 50% probability to be in 
the positive or negative direction. The rates at which 
elementary moves are attempted is thus equal for inte- 
rior monomers and those at the ends. If the first or last 
monomer is selected, similar statements suffice to update 
(Ax,Bx) or (Ajv-ij-Bjv-i)) respectively. The following 
statements can be used to update the end- monomer 1, 
trying to displace it in the negative or positive direction 
based on the value of r% . 

A[ = -(BiVn) 

B[ = (-Ax) An W 

However, in case the first monomer is selected, also its 
coordinate x% needs to be updated. Since this coordinate 
is an integer number that can take a wide range of values, 
its update is not implemented in a multispin coding fash- 
ion. Luckily, the probability to select the first monomer 
decreases with increasing polymer length. 

The motivation for multispin coding lies in its effi- 
ciency. The simulation of the dynamics of 64 polymers 
involves 128-/V elementary moves per time unit. With 
multispin coding, it involves that the above sequence of 
logical operations for interior monomers should be per- 
formed statistically N — 2 times, which involves only 11 
logical operations, 4 loads, and 4 stores. On a fast work- 
station, this takes 41 ns of cpu-time, or 0.32 ns per el- 
ementary move. For updating the last monomers, even 
less logical operations suffice. For the first monomers, 
however, the update of the 64 values of x\ cannot be 
achieved as efficiently, and a loop over all 64 polymers is 
inevitable, with the consequence that the simulation of 
the first monomer requires 5.8 ns cpu-time per elemen- 
tary move. 

The multispin implementation outlined above will thus 
perform 64 simulations in parallel. These simulations are 
however correlated, since they share the sequence of se- 
lected monomers. In fact, if at some point two simula- 
tions are in identical polymer configurations, they will 
stay in identical configurations ever after; in long sim- 
ulations of small systems, this "locking" will inevitably 
happen. Complete locking is avoided by an uncorrelated 
choice in the directions in which end monomers attempt 
to move, using a 64-bit random bit pattern rather than a 
binary choice between all-up or all-down. If desired, more 
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TABLE I: The twelve vectors pointing to nearest neighbors 
of a FCC lattice, expressed in combinations of t, u, v and 
to = — (t + u + v). 
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de-locking can be obtained at the expense of a lower effi- 
ciency, by rejecting a fraction of the allowed moves, also 
using a 64-bit random bit pattern. An important remark 
is, however, that as long as each of the 64 simulations 
are correct in themselves, one obtains 64 unbiased re- 
sults; one should just be careful to assign significance to 
the spread in those 64 simulations. 

B. model for polymer melts 

While monomers in the projected repton model live 
on a one-dimensional lattice, the monomers in the model 
that we propose for polymer melts live on a FCC lat- 
tice. It is helpful to note that the three-dimensional 
hyperplane, located in a four-dimensional hypercubic 
space through the origin and with perpendicular vector 
(1,1,1,1), is such an FCC lattice. Stated differently, 
the set of points x — (a,b,c, d) with integer- valued co- 
ordinates, constrained to a + b + c + d = 0, forms a 
FCC lattice. The twelve vectors pointing to nearest- 
neighbor sites are t = (—1,1,0,0), u = (0,-1,1,0), 
V = (0,0,-1,1), w = -(t + u + v) = (1,0,0,-1), and 
some of their combinations, as listed in Table The 
vector (0, 0, 0, 0) is used as the representation for a zero- 
length bond. 

As in the projected repton model, the state of a poly- 
mer on our FCC lattice can be specified by the location 
j?i of the first monomer, plus additionally the direction 
in which adjacent neighbors are located. We choose for 
the polymer melt model to store the four directional bits 
in a single word: bits k, fc+16, k + 32 and fc+48 of the 64- 
bit word Di indicate the vector pointing from monomer 
i to monomer i + 1. Thus, sixteen polymers are updated 
simultaneously. 

The coordinates of the i th monomer in polymer k can 
then be retrieved by summing over all words Di bits 
k, k + 16, k + 32 and k + 48, yielding respectively the 
numbers a, b, c and d; the monomer position is then 



Xi = at + bu + cv + dw. Note that these summations 
require only 3« operations (a right shift of k bits, mask- 
ing the direction bits, and adding to the sum), since the 
summation in the different bits can be done in a single 
operation. Since in the polymer melt model we need the 
coordinates frequently (every time we attempt a possible 
sideways move), and since our polymers are often several 
hundred monomers long, we do not keep track explicitly 
of only the position x\ of the first monomer, but also of 
the last monomer as well a few other monomers along 
the chain, such that the distance along the chain to a 
monomer with known position is always less than 15. 
Of course, to retrieve the position Xi, we start from the 
nearest monomer with known position in either direction 
along the polymer. 

The implementation proceeds analogous to the pro- 
jected repton model. Also here, if an interior monomer 
can move in one direction, its move in the other direc- 
tion is blocked; this can be exploited as in the repton 
model, by combining two elementary moves. The precise 
operations are: 



2/o 


= Di V (Di > 32) 


2/1 


= 2/ V(y >16) 


zo 


= Dj V (Dj > 32) 


Zl 


= z V {z > 16) 
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= (yi © zi) A (2 16 - 1) 


mi 
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= Di®m x 


D 'j 


= D 3 © mi 



(3) 



Here, A <C k stands for shifting the word A over k bits 
to the left, and M = 2° + 2 16 + 2 32 + 2 48 is a constant, 
used to duplicate the low 16 bits in the higher bits of 
the word. Thus, with only 15 operations, 2 loads, and 
2 stores, we have performed 32 elementary moves. On 
a fast workstation, the above implementation requires 
1.25 ns cpu-time per elementary move. 

As in the projected repton model, if the first monomers 
are displaced, whose positions are tracked, these posi- 
tions have to be updated in a loop over the 16 poly- 
mers. We succeeded in implementing moves of the first 
monomers in 82 ns cpu-time per elementary move. Other 
monomers whose position is tracked require roughly the 
same computational effort. Displacement of the first and 
last monomers is attempted at twice the rate of the other 
monomers, for the same reason as in the projected repton 
model. 

Besides the reptation moves, the dynamics consists of 
sideways moves. If a sideways move is attempted on 
monomer i, the relevant bit patterns are those indicating 
the direction from monomer i— 1 to i and from i to i + 1. 
These bit patterns, as listed in Tabled are all numbers in 
the range to 15, except for the values 5, 10, and 15 that 
do not occur; the number denotes stored length, while 
the other 12 numbers denote bonds to the twelve nearest- 
neighbor sites. For every combination (Di,Dj) of those 
bit patterns, we have pre-computed the lists of all bit pat- 
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terns D[ and Dj after a sideways move. Depending on 
the combination (D,, Dj), at most four different sideways 
moves can be proposed. We therefore have precomputed 
four such lists for D[ and for Dj . In one step, we first se- 
lect randomly the monomer number i and the list number 
k; then we attempt a sideways move of monomer i to the 
position determined by the fe th list; finally, if this move 
does not lead to overlapping monomers, it is accepted. 
This check for overlap requires computing the position 
of monomers, which requires computing the distance to 
the nearest tracked monomer. In our implementation, in 
which this distance is at most 15 monomers, the total 
cpu-time required per such step equals 90 ns. 

Not all moves are attempted with the same frequency. 
Since the long-time dynamics is determined by reptation, 
as argued by de Gennes, the time scale is set such that 
reptation moves in the interior of the polymers are at- 
tempted with unit rate. Sideways moves which do not 
increase the amount of stored length are attempted with 
some rate r s ; the most natural choice for this rate would 
be unit rate once more, but because of the much higher 
computational cost for these moves, we often chose some 
value of r s < 1. As discussed above, consistency with the 
density of stored length of the projected repton model de- 
mands that moves in which the amount of stored length 
is increased are attempted with a rate of 2r s /z. Since 
the most mobile monomers hop away to other sites with 
a total rate of 2r s , and since 16 monomers are potentially 
moved in each step, r s /8 steps as described above should 
be performed per monomer and per unit of time. 

Moves in which the first or last monomer joins its 
neighbor along the chain can be viewed partly as rep- 
tation moves, partly as sideways movement. On these 
grounds, we have decided to attempt these moves with 
rate 1 + r s . Consistency with respect to the density of 
stored length requires that moves in which the first or last 
monomer leaves its neighbor along the chain in a specific 
direction are attempted with a rate of 2(1 + r s )/z. 



IV. APPLICATION: PHASE SEPARATION OF 
A BINARY POLYMER MIXTURE 

To illustrate the efficiency of the above computational 
approach, we performed a simulation of the phase sepa- 
ration of a binary polymer mixture with polymer types 
A and B. The A and B polymers interact with a short- 
range repulsion, described by the Hamiltonian 

H = jJ2 A ) 6 (< J P' , B ) + B ) S (^' > A )> ( 4 ) 

where the summation runs over all pairs of nearest- 
neighbor sites, and a(f) is A, B, or if site f is occupied 
by a polymer of type A, B or empty, respectively. The 
repulsion between A and B polymers provides the driving 
force for the phase separation. 

We simulated a system containing in total 46, 080 poly- 
mers of length L = 100 on a lattice of N = 13, 824, 000 




FIG. 3: Two-dimensional snapshot of a three-dimensional 
simulation of a binary polymer mixture with 46, 080 polymers 
of length L = 100 on a lattice of N = 13, 824, 000 sites, at 
inverse temperature /3J — 0.1. Snapshot a) is taken at t = 0, 
i.e., the equilibrated system at infinite temperature taken as 
the initial configuration. Snapshots b) and c) are taken at 
t = 4.7 ■ 10 5 and 3.7 ■ 10 6 MC time units. 



sites, at inverse temperature 0J — 0.1. The system 
evolves in time through reptation moves, at unit rate, 
in combination with sideways moves at a rate of r s = 
1/30. Figure |3 shows two-dimensional slices of the three- 
dimensional system at times t = 0, t = 4.7 • 10 5 , and 
t = 3.7 • 10 6 MC time units. 

At various times, we determine the two-point correla- 
tion function 



9AB{r) 



N 



J2f* S(ap,A) S(ap +? , B) 



(5) 



From this function we determine the spherically-averaged 
radial distribution function, defined as RDF(r) = 1 — 
g(r). This function is 1 at r = 0, then decreases, and 
eventually approaches for large r. After some time, 
the conserved dynamics gives rise to damped oscillations 
in the RDF. The frequency of these oscillations can be 
determined from the shortest distance ro at which the 
RDF equals zero. The typical domain size d(t) is then 
obtained as twice this distance. 
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Cubic root of simulation time (x 1000 MC time units) 



FIG. 4: Domain size of a phase-separating binary polymer 
mixture as function of the cubic root of the simulation time. 
The three symbols indicate three independent simulations. A 
good agreement with domain growth d(t) ~ t 1 ^ 3 is found. 

Figure^Jshows the domain size of the phase-separating 
mixture at different times. The domain size d(t), as de- 
fined above, is plotted as a function of the cubic root 
of the simulation time in Monte Carlo sweeps. The 
straight line shows that the domain size grows with 



d(t) ~ t 1 / 3 . This is to be expected in a system with 
conserved order dynamics and without hydrodynamics, 
as in "Model B" [2f|. 

V. SUMMARY AND OUTLOOK 

In summary, we have presented a lattice model to 
study the dynamics of polymeric systems. Compared to 
other lattice polymer models in current use, this model 
lends itself for highly efficient simulation: employing mul- 
tispin coding, the computer time required per elementary 
move is reduced by more than three orders of magnitude, 
to a few nanoseconds. 

To demonstrate the strength of our model, we simu- 
lated the phase separation of a binary mixture of poly- 
mers. We performed simulations of a system with about 
50, 000 polymers, each containing 100 monomers, located 
on a lattice with about 14 million sites. We verified that 
over the final two decades in time, the domain size d(t) 
grows according to d(t) ~ i 1 / 3 , as expected for a system 
with overdamped dynamics and a local conservation law. 
This simulation involves 3.4 • 10 13 elementary moves, but 
could nevertheless be carried out on a single-processor 
workstation in about 20 days. 

Recently, we have used this lattice polymer model to 
study fractionation in quasi-binary (polydisperse) poly- 
mer mixtures |27^ . Currently, we are using the model 
presented here to study the dynamics of absorption of a 
polydisperse polymer mixture on a surface, as well as the 
sieving process of polydisperse polymer mixtures through 
nanopores. 
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